The Signature of the Ice Line and Modest Type I Migration in 
the Observed Exoplanet Mass-Semimajor Axis Distribution 

Kevin C. Schlaufman^ and D. N. C. Lin^ 

Astronomy and Astrophysics Department, University of California, Santa Cruz, CA 95064; 
. . . kcs@ucolick.org and lin@ucolick.org 

00 ; 

o : 

^ ■ S. Ida 

^ ■ Tokyo Institute of Technology, Ookayama, Meguro-ku, Tokyo 152-8551, Japan; 
■ ida@geo.titech.ac.jp 

ABSTRACT 

2 I Existing exoplanet radial velocity surveys are complete in the planetary mass- 

^ \ semimajor axis [Mp — a) plane over the range 0.1 AU < a < 2.0 AU where 

li^l Mp > 100 M^. We marginalize over mass in this complete domain of parame- 

I I ter space and demonstrate that the observed a distribution is inconsistent with 

I models of planet formation that use the full Type I migration rate derived from a 

if^ \ linear theory and that do not include the effect of the ice line on the disk surface 

^ I density profile. However, the efficiency of Type I migration can be suppressed 

^ I by both nonlinear feedback and the barriers introduced by local maxima in the 

O ' disk pressure distribution, and we confirm that the synthesized M„ — a distri- 
00 ' 

^ I bution is compatible with the observed data if we account for both retention of 

^ I protoplanetary embryos near the ice line and an order-of-magnitude reduction 

I in the efficiency of Type I migration. The validity of these assumption can be 

H ' checked because they also predict a population of short-period rocky planets with 



a range of masses comparable to that of the Earth as well as a "desert" in the 
Mp — a distribution centered around Mp ~ 30 — 50 M® and a < 1 AU. We show 
that the expected "desert" in the Mp — a plane will be discernible by a radial 
velocity survey with 1 m s~^ precision and n ~ 700 radial velocity observations 
of program stars. 

Subject headings: planetary systems — planetary systems: formation — plane- 
tary systems: protoplanetary disks 
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Introduction 



Over 200 planets with reliable mass (Mp) and semimajor axis (a) measurements have 
been discovered around nearby FGK stars in the past decade. At the same time, attempts to 
build a comprehensive deterministic theory of planet formation have lead to the de velopment 
of pop ulation synthesis models based on the sequential accretion scenario. In llda &: Lin 
( I2OO4J ). two of us studied the growth of planetesimals into dynamically isolated embryos 
as well as their tidal interactions with their parent disks. Using the observed ranges of 
disk mass, size, and accretion rate we showed that a fraction of embryos evolve into cores 
with more than a few Earth masses (M^), accrete massive envelopes, open up gaps near 
their orbits, and attain asymptotic masses comparable to that of Jupiter. In some massive 
and persistent disks, the newly formed gas giant planets may migrate toward the proximity 
of their host stars. In the end these simulations produced the distribution of dynamical 
and structural properties of planets. Presently, the observed sample of extrasolar planetary 
properties has become large enough to enable direct comparisons between the theoretically 
predicted and observed Mp — a distributions that not only delineate the dominant physical 
mechanisms at work in planet formation, but also provide quantitative constraints on the 
efficiencies of those processes. 

In the latest update of the planet forma tion models two of us have incorporated the 
effect of Type I migration (llda fc LinI l2008al ). This process is a direct consequence of a 
protoplanetary core's tidal interaction with its parent protopl anetary disk. The efficiency 



of Type I migration was first d etermined by a linear theory (iGoldreich fc Tremaind Il980 



WardI Il986l : iTanaka et al.l |2002| ) that neglected the embryo's perturbation on the surface 
density distribution of its parent protoplanetary disk. In the environment of a minimum 
mass nebula though, this efficiency factor would imply that a protoplanetary embryo with 
mass a fraction of an Earth mass would migrate from ~ 1 AU into its host star before the 
severe depletion of the disk gas that is known to occur over a time scale of several Myr. 
Although this critical mass at which Type I migration causes AU-scale migration on the disk 
depletion time increases with a, it is still difficult to retain sufficiently massive cores for the 
onset of dynamical acc retion of gas. This argument implies that gas giants should be very 
rare (llda fc Linll2008al ) and this paradox has led to many in-depth analyses of the Type I 
migration process. Numerica l nonlinear simulations of Type I migration were reviewed by 
Papaloizou fc TerquemI (120061 ) and many potentia l explanations for slow Type I migration are 
in the literature: intrinsic turbulence in the disk (iLaughlin et al.ll2004l : iNelson fc Papaloizou 



20041 ). self-induced unst able flow (IKo 



hydrodynamic feedback (iMasset et a 



ture gradients (IMasset et al.ll2006bl ). 



er fc 



2004 



Li et al.l 20051). nonlinear radiative and 



200631). and variat ion in surface density and tempera- 



Dobbs-Dixonl (120071 ) has shown that in some situations 



the nonlinear Type I migration rate can be less than 10% of the linear prediction. 
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Another issue studied by two of us (llda &: Linll2008bl . IL hereafter) is the critical embryo 
mass Merit > at least a few required by current models for runaway gas accretion. The 
embryo is limited by its isolation mass Mj^o, and the solid surface density profile of the 
minimum mass solar nebula (MMSN) oc a~^/^ requires that the Mj^o scales like a^^^. 
On the other hand, the timescale for growth Tc^acc scales with a^^/^°. As a result, after 
the characteristic gas depletion time Tdep the most massive embryos near the ice line have 
masses Mc ~ Mj^o < Merit- However, since we observe many exoplanets with Jupiter 
masses a t a ~ 1 AU, there m ust be some physical process that is neglected in this simple 
analysis. iKretke &: LinI (120071 ) outlined one possible solution to this problem in which solids 
are trapped near regions of the disk where the local pressure requires the gas to rotate with 
super-Keplarian velocities. When th e combined coritribut ion of both Lindblad and corotation 
resonances are taken into account (IMasset et al.l l2006bl ) . the migration of protoplanetary 
cores may be suppressed as well (IILI). 



In this paper, we utilize the observed data to calibrate the population synthesis models. 
In §2 we quantitatively show that the existing distribution of exoplanets cannot be explained 
by models of planet formation that apply the full Type I migration rate predicted from linear 
theory and that do not include the effects of the ice line. We also point out that the existing 
observed and synthesized Mp — a distributions are in agreement with each other if we take 
into account the effect of an ice line barrier and assume a reduction in the magnitude of Type 
I migration. In addition, we describe the parameters of a radial velocity survey capable of 
verifying the existence of "desert" in the Mp — a diagram predicted by |IL|. In §3 we consider 
the implications of these models and suggest methods to test our assumptions. In §4 we 
summarize our findings. 



2. Analysis 



We u s e the ideas presented in iNarayan et al.l (120051 ) an d the formalism d evelop ed in 



Cummingl (12004 ) to approximately reproduce the result of ICumming et al.l (120081 ) that 



showed that the current California and Carnegie Planet Search (CCPS) has announced 
all planets with orbital period P < 2000 days, stellar reflex velocities K > 20 m s~^, and 
eccentricities e < 0.6. In particular, we determine which of the simulated exoplanet systems 
from |llj would be detectable by a radial velocity survey with precision and cadence similar 
to present-day radial velocity surveys like the CCPS and the High Accuracy Radial Velocity 
Planetary Search Project (HARPS). In this approach, we combine Equation (26) and (30) of 
Cummingl (120041 ) and declare that all the synthesized planetary systems from |lLj with mass 



Mp > M50 where 
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are detectable. In Equation ([T]), is the number of radial velocity observations, a is the 
RMS of spectrograph precision and stellar jitter, i is the inclination of the exoplanet's orbit, 
P is the orbital period in years, and M^, is the host stellar mass in solar masses. In the limit 
of large N, M ^ 100 is the number of independent frequencies searched and F ^ 0.01 is 
the false alarm probability - the numerical values are correct to order-of-magnitude and in 
any case they only very weakly influence our estimate of M50. We note that E quation (fTI) 
i s form ally correct only for single planetary systems in circular orbits; however, ICumming 
( 120041 ) shows that in the limit of large N Equation ([1]) applies to multiple planet systems 
and in the case where e < 0.6. Therefore, we set M50 = cx) if e > 0.6. We assume a velocity 
resolution of 1 m s~^. 

We then carry out a Monte Carlo simulation in which we assign each of the simulated 
planetary systems from |llj a random host stellar mass, stellar jitter, eccentricity, inclina- 
tion, and number of radial velocity observations. In this prescription, we use the empirical 
distributions for host stellar mass, stellar jitter , eccen tricity, and number of radial velocity 
observations given in the updated Butler et al. (2006) catalo jl| of all known exoplanets; we 
use the standard distribution for random inclinations. We compute which planetary systems 
have Mp > M50 and declare that these systems are detectable in this iteration. 

We repeat this process 1000 times. We consider all planetary systems that are detectable 
according to the M50 criterion in 90% of the Monte Carlo iterations robustly detectable. 
Averaged over host stellar mass, stellar jitter, eccentricity, inclination, and number of radial 
velocity observations we find that all of the simulated planets from|lLj in the range 0.1 AU 
< a < 2.0 AU with 



Mp > 



174 ( — ) + 47 

au 



(2) 



are robustly detectable. We include the results of this calculation in Figure [H The robustly 
detectable systems are marked with solid circles and we plot the observed Mp — a distribution 
of known planetary systems as the solid squares. 



Maintained at |http: / /www. exoplanets. org 
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ILl generated a set of 12 realizations of the Mp — a plane under different physical assump- 
tions described in Table [H One group of models result from disks with the characteristic 
bump in gas surface density due to the couphng effect of the MRI activity and the ice line 
and with an enhancement in solid surface density E^, another group of models has the bump 
in Eg but not in E^, and the last group ignores the effects the ice line would have on E^ 
and Erf. For each group of models, four different parameterizations of the Type I migration 
rate were used: 100%, 30%, 10%, and 3% of the full rate from linear theory. For each group 
of models we marginalize the 2D Mp — a distribution of the simulated exoplanet systems 
in the complete region over planetary mass, leaving us with ID distributions in semimajor 
axis. We plot histograms for each group of models: Figure [2] corresponds to the models that 
disregard the effects of the ice line on E^ and E^, Figure [3] corresponds to the models that 
includes the effect of the ice line on Eg but not E^, and Figure H] corresponds to the models 
that include the effects of the ice line on Eg and E^. 

We determine the model that best matches the observed Mp — a distribution in the 
complete region by computing for each model the p-value distribution that results from 
comparing 1000 bootstrap resamplings from that model in the complete region with 1000 
resamplings of the observed data in the complete region. We include the results of this 
calculation in Figure [5] and we also report median p-values and 95% intervals in the last 
three columns of Table [TJ We find that only models that include the effects of the ice line 
on both the gas and solid surface density of the disk and apply a Type I migration rate an 
order-of-magnitude less than that predicted by linear theory (models COIC and C003C) are 
consistent with the observed data in the complete region. Models which neglect the presence 
of the ice-line barrier generally do not yield the observed up-turn in the period distribution 
of the known planets and are rejected at very high confidence. Models with efficient Type I 
migration generally under predict the fraction of stars with detectable gas giants, especially 
for those with a outside the ice line. Therefore, we argue that the population synthesis 



models presented in incorporating a Type I migration rate much reduced from linear 
theory and the effects of the ice line on both the solid and gas surface densities (E^ and Eg 
respectively) are at least plausible deterministic models for giant planet formation. We also 
note that the observed data suggests that the ad hoc prescription for the location of the ice 



line used by |IL| underestimates its radius by perhaps even a factor of two. 



Furthermore, in the presentation of their population synthesis models |IL| illustrated the 
presence of a "desert" in the Mp — a distribution. This sparsely populated region is depleted 
due to both Type I migration and runaway gas accretion. For Models COIC and C003C, 
this "desert" is centered around Mp ~ 30 and a < 1 AU. We use the same detection 
strategy described above, only now we model the number of radial velocity observations of 
each planetary system over a period of about ten years as a Gaussian random variable with 
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mean and standard deviation cr„; we round each random deviate to the nearest whole 
number. We then fit a two component Gaussian mixture model to the Mp — a distribution 
of all robustly detected planet with mass Mp < 100 M^, and we say the "desert" is detected 
if the mean vectors of the two Gaussians are offset by more than 0.6 in log a and the minor 
axis of the Gaussian at smaller orbital radius is larger than the minor axis of the Gaussian 
at larger orbital radius. In other words, if the two components of the mixture model bracket 
the corner of a region devoid of extrasolar planets - a metaphorical "desert" - we say the 
"desert" is resolved. We find that when /x„ = 700 and ct„ = 50, the two components of the 
mixture model bracket a barren region and therefore the "desert" is resolved more than 90% 
of the time. We illustrate results of our calculation in Figure [61 As a result, a radial velocity 
campaign with the parameters described above will be able to confirm the fidelity of the 
models presented in 



IL. 



Discussion 



The key prediction of |IL| is that for masses at which the dominant migration mechanism 
is Type I migration - Mp < 50 - there will be a dearth of exoplanets within 1 AU of 
their host stars, simply because the timescale for Type I migration is so much shorter than 
the disk dispersal time. Since they are formed interior to the ice line, these planets are likely 
to be rocky and have mass a few M^. also predict an overdensity of gas giant planets at 



~ 2 AU resulting from the ice line. We show in Figure [6] that one can quantitatively detect 
these features in the observed Mp — a with a radial velocity survey with 1 m s~^ precision and 
about 700 radial velocity observations, or about ten years worth of data from the Automated 
Planet Finder (APF). In future work, we will examine the ability of missions like the Space 
Interferometry Mission (SIM) to verify the same features in Mp — a plane. 

In addition to the upper and lower Mp bound, the "desert" is also surrounded by pop- 
ulated domains in the a distribution. While the ice line provides a strong up-turn at a few 
AU, the models by [llj also imply a large population of short-period rocky planets as a con- 
sequence of Type I migration. Despite an order-of-magnitude decrement in the efficiency of 
Type I migration, the simulated Mp — a distribution of models COIC and C003C indicate 
that in the proximity of their host stars, rocky planets with Mp ~ a few are at least 



an order-of-magnitude more common than close-in gas giants (see Figure 6 of IlLI). Other 
authors have already pointed out the observational difficulties inhe rent in the search fo r this 
population of "super-Earths" through radial velocity observations (INarayan et al.ll2005l ). Fi- 
nally, we note that the models of IL did not include dynamical interactions between planets 
in multiple planetary systems, and these interactions can broaden the simulated a distribu- 
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tion and eccentricity distribution. We will include these effects in future generations of the 
population synthesis models. 



4. Conclusion 

We used the fact that existing exoplanet radial velocity surveys are complete in the 
planetary mass-semimajor axis (Mp — a) plane where 0.1 AU < a < 2.0 AU and Mp is in the 
range specified by Equation ([2]) to show that the observed semimajor axis distribution in the 
complete region cannot be explained by models of planet formation that use the full Type I 
migration rate predicted by linear theory and that do not include the effects of the ice line. 
Moreover, we also demonstrated that the expected "desert" in the Mp — a plane at about 



Mp ~ 30 Mq) and a < 1 AU predicted by |lLj will be discernible by a radial velocity survey 



with 1 m s~^ precision and n ~ 700 radial velocity observations of program stars. Such an 
observational campaign will also verify the predicted inner boundary of the "desert" where 
we expect a large population of super-Earths have migrated to and halted in the proximity 
of their host stars. 

We thank A. Gumming, G. Laughlin, G. Marcy, and Michel Mayor for useful conver- 
sation and the anonymous referee for some insightful comments. This research has made 
use of NASA's Astrophysics Data System Bibliographic Services. This material is based 
upon work supported under a National Science Foundation Graduate Research Fellowship, 
NASA (NAGS5-11779, NNG06-GF45G, NNX07A-L13G, NNX07AI88G), JPL (1270927), 
NSF(AST-0507424), and JSPS. 



REFERENCES 

Butler, R. R, et al. 2006, ApJ, 646, 505 
Gumming, A. 2004, MNRAS, 354, 1165 

Gumming, A., Butler, R. R, Marcy, G. W., Vogt, S. S., Wright, J. T., & Fischer, D. A. 2008, 
PASP, 120, 531 

Dobbs-Dixon, I. M. 2007, Ph.D. Thesis 

Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 

Ida, S., & Lin, D. N. G. 2004, ApJ, 604, 388 



- 8 - 

Ida, S., & Lin, D. N. C. 2008, ApJ, 673, 487 

Ida, S., & Lin, D. N. C. 2008, ArXiv e-prints, 802. larXiv:0802.1114l 
Laughlin, G., Steinacker, A., & Adams, F. C. 2004, ApJ, 608, 489 

Li, H., Li, S., Koller, J., Wendroff, B. B., Liska, R., Orban, C. M., Liang, E. P. T., & Lin, 
D. N. C. 2005, ApJ, 624, 1003 

Koller, J., & Li, H. 2004, The Search for Other Worlds, 713, 63 

Kretke, K. A., & Lin, D. N. C. 2007, ApJ, 664, L55 

Masset, F. S., D'Angelo, G., & Kley, W. 2006, ApJ, 652, 730 

Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 642, 478 

Narayan, R., Gumming, A., & Lin, D. N. G. 2005, ApJ, 620, 1002 

Nelson, R. R, & Papaloizou, J. G. B. 2004, MNRAS, 350, 849 

Papaloizou, J. G. B., & Terquem, G. 2006, Reports of Progress in Physics, 69, 119 

Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 

Ward, W. R. 1986, Icarus, 67, 164 



This preprint was prepared with the AAS IfTJrjX macros v5.2. 



- 9 - 



Table 1. Model Descriptions from llda fc LinI (j2008bl ) 



Name 




Hg Enhanced 


Srf Enhanced 






pC 








CIC 


1 


Yes 


Yes 


6.9 X 10- 


-6 


2.8 X 10- 


-3 


1.7 X 10- 


-1 


C03C 


0.3 


Yes 


Yes 


1.6 X 10- 


-6 


1.7 X 10- 


-3 


1.2 X 10- 


-1 


COIC 


0.1 


Yes 


Yes 


9.9 X 10- 


-3 


8.6 X 10- 


-2 


7.2 X 10- 


-1 


C003C 


0.03 


Yes 


Yes 


3.2 X 10- 


-3 


6.0 X 10- 


-2 


6.0 X 10- 


-1 


CIB 


1 


Yes 


No 


2.9 X 10- 


12 


7.3 X 10- 


-8 


9.6 X 10- 


-5 


C03B 


0.3 


Yes 


No 


5.8 X 10- 


-9 


1.4 X 10- 


~5 


2.8 X 10- 


-3 


COIB 


0.1 


Yes 


No 


3.4 X 10- 


-6 


2.8 X 10- 


-3 


2.2 X 10- 


-1 


C003B 


0.03 


Yes 


No 


5.1 X 10- 


-5 


1.2 X 10- 


-2 


3.8 X 10- 


-1 


Cl_p4 


1 


No 


No 


1.0 X 10- 


16 


1.0 X 10- 


16 


1.0 X 10- 


-16 


C03_p4 


0.3 


No 


No 


1.0 X 10- 


16 


1.0 X 10- 


16 


1.0 X 10- 


-16 


C01_p4 


0.1 


No 


No 


1.0 X 10" 


16 


1.0 X 10- 


16 


1.0 X 10- 


-16 


C003_p4 


0.03 


No 


No 


4.4 X 10" 


16 


2.2 X 10- 


11 


3.6 X 10- 


-7 



^From llda fc LinI (j2008bl ) - Ci equivalent to the fraction of the full Type I migration 
rate predicted from linear theory applied during the simulation 

'^ower bound of an interval centered on median p-value which contains 95% of our boot- 
strap realizations 

■^Median p-value of our bootstrap realizations 

"^Upper bound of an interval centered on median j>value which contains 95% of our 
bootstrap realizations 
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Fig. 1. — Results of our Mo n te Carl o simulation. The open circles are the simulated planets 
from Figure 3 of llda fc Liru (j2008bl ): the model includes the effects of the snow line on the 
surface density of gas and dust T,d- The Type I migration rate in the simulation is 
0.1 times the prediction from linear theory. The filled circles are planets that are robustly 
detected, that is, planets that would be detected at least 90% of the time by current radial 
velocity surveys. We also plot all known exoplanet planets as filled squares. All simulated 
planets in the range 0.1 AU < a < 2.0 AU with Mp as specified by Equation ([2]) denoted by 
the heavy black lines are robustly detected, so existing radial velocity surveys are complete 
in that range. See the electronic edition of the Journal for a color version of this figure. 
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Fig. 2. — ID distribution in semimajor axis derived from the p r ojectio n of the 2D distribution 
in the complete region of the Mp — a plane from llda &: LinI (j2008bl ) for the disks without 
the bump in T,g or S^. The upper left panel uses the full Type I migration rate from linear 
theory, the upper right panel uses 30% of the full rate, the bottom left panel uses 10% of 
the full rate, and the bottom right panel uses 3% of the full rate. See the electronic edition 
of the Journal for a color version of this figure. 
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Simulation C1B 
Observation 



H 1 1 I I I I 



Simulation C03B 
Observation 



H 1 — I I I I I 



Simulation C01B 
Observation 



Simulation C003B 
Observation 



10" 
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10" 



10^ 
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Fig. 3. — ID distribution in semimajor axis de rived from the pro j ection of the 2D distribution 
in the complete region of the Mp—a plane from llda fc Liru (j2008bl ) for the disks with the bump 
in Eg but without enhancement in E^. The upper left panel uses the full Type I migration 
rate from linear theory, the upper right panel uses 30% of the full rate, the bottom left panel 
uses 10% of the full rate, and the bottom right panel uses 3% of the full rate. The error bars 
indicate the 2-a region. See the electronic edition of the Journal for a color version of this 
figure. 
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Simulation C1C 
Observation 



Simulation C03C 
Observation 
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Simulation C01C 
Observation 



Simulation C003C 
Observation 



10" 



a[AU] 



10" 



10^ 



a[AU] 



Fig. 4. — ID distribution in semimajor axis derive d from the projectio n of the 2D distribution 
in the complete region of the Mp — a plane from llda fc LinI (j2008bl ) for the disks with the 
bump in both Eg and S^. The upper left panel uses the full Type I migration rate from 
linear theory, the upper right panel uses 30% of the full rate, the bottom left panel uses 10% 
of the full rate, and the bottom right panel uses 3% of the full rate. The error bars indicate 
the 2-0" region. See the electronic edition of the Journal for a color version of this figure. 
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Fig. 5. — Distribution of the Kolmogorov-Smirnov test p-values resulting from 1000 boot- 
strap resamplings; similar distributions will have a sharply-peaked p- value distributions with 
a maximum near p ~ 1 ^ logp ~ 0. The top row shows the results for models including 
the bump in both Eg and S^, the middle row shows models including just the bump in Eg, 
and the bottom row shows models with no bump in Eg or E^. For all rows, the first column 
shows models with the full Type I migration rate, the second column shows models with 30% 
of the full Type I migration rate, the third column shows models with 10% of the full Type 
I migration rate, and the fourth column shows models with 3% of the full Type I migration 
rate. The p- value distribution of the model with the bump in both Eg and E^ and 10% of the 
full Type I migration rate is the best match to the observed data. There are no histograms 
for models Cl_p4, C03_p4, or C01_p4 because the p-values were vanishingly small. See the 
electronic edition of the Journal for a color version of this figure. 
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Fig. 6. — Results of our Monte Carlo simulation. The open circles are the simulated planets 
of model COIC. As in Figure [H the filled circles are planets that are robustly detected by 
a radial velocity survey with n ~ 700 observations per program star; we plot all known 
exoplanet planets as filled squares. The centers of the best fit Gaussian mixture to all 
robustly detected exoplanets below the dashed line at 100 models are marked by the 
two X's, while the characteristic ellipses of each component are the solid lines. We say the 
"desert" is detected if the mean vectors of the two components are offset by more than 0.6 
in log a and the minor axis of the ellipse at small orbital radius is larger than the minor axis 
of the ellipse at large orbital radius. We find that the "desert" is detected more than 90% 
of the time when = 700. See the electronic edition of the Journal for a color version of 
this figure. 



